// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_SPARSEPRODUCT_H
#define EIGEN_SPARSEPRODUCT_H

namespace Eigen {

/** \returns an expression of the product of two sparse matrices.
  * By default a conservative product preserving the symbolic non zeros is performed.
  * The automatic pruning of the small values can be achieved by calling the pruned() function
  * in which case a totally different product algorithm is employed:
  * \code
  * C = (A*B).pruned();             // suppress numerical zeros (exact)
  * C = (A*B).pruned(ref);
  * C = (A*B).pruned(ref,epsilon);
  * \endcode
  * where \c ref is a meaningful non zero reference value.
  * */
template <typename Derived>
template <typename OtherDerived>
inline const Product<Derived, OtherDerived, AliasFreeProduct> SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived>& other) const
{
    return Product<Derived, OtherDerived, AliasFreeProduct>(derived(), other.derived());
}

namespace internal {

    // sparse * sparse
    template <typename Lhs, typename Rhs, int ProductType> struct generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType>
    {
        template <typename Dest> static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
        {
            evalTo(dst, lhs, rhs, typename evaluator_traits<Dest>::Shape());
        }

        // dense += sparse * sparse
        template <typename Dest, typename ActualLhs>
        static void addTo(Dest& dst,
                          const ActualLhs& lhs,
                          const Rhs& rhs,
                          typename enable_if<is_same<typename evaluator_traits<Dest>::Shape, DenseShape>::value, int*>::type* = 0)
        {
            typedef typename nested_eval<ActualLhs, Dynamic>::type LhsNested;
            typedef typename nested_eval<Rhs, Dynamic>::type RhsNested;
            LhsNested lhsNested(lhs);
            RhsNested rhsNested(rhs);
            internal::sparse_sparse_to_dense_product_selector<typename remove_all<LhsNested>::type, typename remove_all<RhsNested>::type, Dest>::run(
                lhsNested, rhsNested, dst);
        }

        // dense -= sparse * sparse
        template <typename Dest>
        static void subTo(Dest& dst,
                          const Lhs& lhs,
                          const Rhs& rhs,
                          typename enable_if<is_same<typename evaluator_traits<Dest>::Shape, DenseShape>::value, int*>::type* = 0)
        {
            addTo(dst, -lhs, rhs);
        }

    protected:
        // sparse = sparse * sparse
        template <typename Dest> static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, SparseShape)
        {
            typedef typename nested_eval<Lhs, Dynamic>::type LhsNested;
            typedef typename nested_eval<Rhs, Dynamic>::type RhsNested;
            LhsNested lhsNested(lhs);
            RhsNested rhsNested(rhs);
            internal::conservative_sparse_sparse_product_selector<typename remove_all<LhsNested>::type, typename remove_all<RhsNested>::type, Dest>::run(
                lhsNested, rhsNested, dst);
        }

        // dense = sparse * sparse
        template <typename Dest> static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, DenseShape)
        {
            dst.setZero();
            addTo(dst, lhs, rhs);
        }
    };

    // sparse * sparse-triangular
    template <typename Lhs, typename Rhs, int ProductType>
    struct generic_product_impl<Lhs, Rhs, SparseShape, SparseTriangularShape, ProductType>
        : public generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType>
    {
    };

    // sparse-triangular * sparse
    template <typename Lhs, typename Rhs, int ProductType>
    struct generic_product_impl<Lhs, Rhs, SparseTriangularShape, SparseShape, ProductType>
        : public generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType>
    {
    };

    // dense = sparse-product (can be sparse*sparse, sparse*perm, etc.)
    template <typename DstXprType, typename Lhs, typename Rhs>
    struct Assignment<DstXprType,
                      Product<Lhs, Rhs, AliasFreeProduct>,
                      internal::assign_op<typename DstXprType::Scalar, typename Product<Lhs, Rhs, AliasFreeProduct>::Scalar>,
                      Sparse2Dense>
    {
        typedef Product<Lhs, Rhs, AliasFreeProduct> SrcXprType;
        static void run(DstXprType& dst, const SrcXprType& src, const internal::assign_op<typename DstXprType::Scalar, typename SrcXprType::Scalar>&)
        {
            Index dstRows = src.rows();
            Index dstCols = src.cols();
            if ((dst.rows() != dstRows) || (dst.cols() != dstCols))
                dst.resize(dstRows, dstCols);

            generic_product_impl<Lhs, Rhs>::evalTo(dst, src.lhs(), src.rhs());
        }
    };

    // dense += sparse-product (can be sparse*sparse, sparse*perm, etc.)
    template <typename DstXprType, typename Lhs, typename Rhs>
    struct Assignment<DstXprType,
                      Product<Lhs, Rhs, AliasFreeProduct>,
                      internal::add_assign_op<typename DstXprType::Scalar, typename Product<Lhs, Rhs, AliasFreeProduct>::Scalar>,
                      Sparse2Dense>
    {
        typedef Product<Lhs, Rhs, AliasFreeProduct> SrcXprType;
        static void run(DstXprType& dst, const SrcXprType& src, const internal::add_assign_op<typename DstXprType::Scalar, typename SrcXprType::Scalar>&)
        {
            generic_product_impl<Lhs, Rhs>::addTo(dst, src.lhs(), src.rhs());
        }
    };

    // dense -= sparse-product (can be sparse*sparse, sparse*perm, etc.)
    template <typename DstXprType, typename Lhs, typename Rhs>
    struct Assignment<DstXprType,
                      Product<Lhs, Rhs, AliasFreeProduct>,
                      internal::sub_assign_op<typename DstXprType::Scalar, typename Product<Lhs, Rhs, AliasFreeProduct>::Scalar>,
                      Sparse2Dense>
    {
        typedef Product<Lhs, Rhs, AliasFreeProduct> SrcXprType;
        static void run(DstXprType& dst, const SrcXprType& src, const internal::sub_assign_op<typename DstXprType::Scalar, typename SrcXprType::Scalar>&)
        {
            generic_product_impl<Lhs, Rhs>::subTo(dst, src.lhs(), src.rhs());
        }
    };

    template <typename Lhs, typename Rhs, int Options>
    struct unary_evaluator<SparseView<Product<Lhs, Rhs, Options>>, IteratorBased> : public evaluator<typename Product<Lhs, Rhs, DefaultProduct>::PlainObject>
    {
        typedef SparseView<Product<Lhs, Rhs, Options>> XprType;
        typedef typename XprType::PlainObject PlainObject;
        typedef evaluator<PlainObject> Base;

        explicit unary_evaluator(const XprType& xpr) : m_result(xpr.rows(), xpr.cols())
        {
            using std::abs;
            ::new (static_cast<Base*>(this)) Base(m_result);
            typedef typename nested_eval<Lhs, Dynamic>::type LhsNested;
            typedef typename nested_eval<Rhs, Dynamic>::type RhsNested;
            LhsNested lhsNested(xpr.nestedExpression().lhs());
            RhsNested rhsNested(xpr.nestedExpression().rhs());

            internal::sparse_sparse_product_with_pruning_selector<typename remove_all<LhsNested>::type, typename remove_all<RhsNested>::type, PlainObject>::run(
                lhsNested, rhsNested, m_result, abs(xpr.reference()) * xpr.epsilon());
        }

    protected:
        PlainObject m_result;
    };

}  // end namespace internal

// sparse matrix = sparse-product (can be sparse*sparse, sparse*perm, etc.)
template <typename Scalar, int _Options, typename _StorageIndex>
template <typename Lhs, typename Rhs>
SparseMatrix<Scalar, _Options, _StorageIndex>& SparseMatrix<Scalar, _Options, _StorageIndex>::operator=(const Product<Lhs, Rhs, AliasFreeProduct>& src)
{
    // std::cout << "in Assignment : " << DstOptions << "\n";
    SparseMatrix dst(src.rows(), src.cols());
    internal::generic_product_impl<Lhs, Rhs>::evalTo(dst, src.lhs(), src.rhs());
    this->swap(dst);
    return *this;
}

}  // end namespace Eigen

#endif  // EIGEN_SPARSEPRODUCT_H
